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1  Objectives 


Reacting  flow  simulations  can  be  significantly  improved  through  the  utilization  of  detailed  chem¬ 
istry  mechanisms.  Two  classes  of  problems  that  may  benefit  include  chemical  laser  systems  and 
hydrocarbon-fueled  propulsion  systems;  each  of  which  typically  involve  on  the  order  of  50  to  200 
chemical  species,  200  rate  equations,  and  time  scales  ranging  from  10~9  seconds  to  over  1  second. 
In  general,  implicit  CFD  algorithms  are  necessary  to  bypass  chemical  stiffness,  but  these  schemes 
become  very  computationally  demanding  (both  in  terms  of  required  memory  and  operation  count) 
as  the  number  of  species  grows.  As  a  result,  incorporation  of  detailed  chemistry  mechanisms  into 
geometrically  complex,  turbulent  CFD  simulations  becomes  computationally  prohibitive:  Current 
approaches  rely  on  reduced  chemical  mechanisms  ( e.g 20  to  30  species  or  less)  which  require  con¬ 
siderable  time  and  effort  to  develop  and  become  less  accurate  with  each  simplification. 

Our  approach  employs  a  loosely-coupled  fluid-dynamic/chemistry  formulation,  that  utilizes  implicit 
solution  methods  for  the  fluid-dynamic  equations  and  low-cost  explicit  or  semi-implicit  diagonal 
methods  for  the  species  continuity  equations.  The  solution  scheme  for  the  chemical  system  ( i.e 
the  species  continuity  equations)  utilizes  a  time-operator  splitting  method  and  a  variable  coefficient 
ordinary  differential  equation  solver  (DVODE)  for  the  reaction  fractional  step  to  appropriately  handle 
chemical  stiffness.  The  DVODE  solver  appears  to  be  very  efficient,  and  this  approach  allows  the  use 
of  less-expensive  explicit  or  semi-implicit  solution  algorithms  for  the  species  transport  fractional  step. 
Our  ultimate  goal  is  to  generalize,  implement,  and  demonstrate  (using  our  GASP  CFD  flow  solver) 
that  this  scheme  can  be  used  to  significantly  reduce  the  costs  associated  with  simulating  complex 
flows  with  detailed  chemical  mechanisms. 

2  Status  of  Effort 

To  date,  we  have  developed  and  implemented  (into  our  GASP  CFD  flow  solver)  a  loosely-coupled 
strategy  for  the  fluid-dynamic  and  chemical  systems.  The  approach  uses  the  lower-upper,  symmetric, 
Gauss-Siedel  (LU-SGS)  implicit  scheme  to  advance  the  fluid-dynamic  system.  In  most  cases,  the 
chemical  reactions  have  much  smaller  time  scales  than  those  associated  with  the  flow,  resulting  in 
stiffness  due  to  the  source  terms.  Solution  of  the  chemical  system  uses  a  time-operator  splitting 
approach  that  allows  isolation  of  the  stiff  source  term  from  the  transport  terms.  As  a  result,  the 
solution  process  can  be  divided  into  two  fractional  steps  -  a  transport  fractional  step  and  chemical 
production  (or  reaction)  fractional  step.  The  chemical  transport  fractional  step  is  currently  solved 
using  the  simple,  low-cost  Euler  explicit  scheme.  The  chemical  production  fractional  step  involves 
solving  a  system  of  ordinary  differential  equations  representing  the  chemical  source  term  and  is 
carried  out  using  a  variable  coefficient  ordinary  differential  equation  solver  (DVODE).  The  DVODE 
solver  has  been  integrated  into  the  GASP  CFD  flow  solver  and  appears  to  be  very  efficient.  We 
have  studied  a  number  of  steady-state  flow  problems,  examining  the  usefulness  of  this  scheme  for 
reducing  chemical  stiffness  and  lowering  computational  costs.  Our  results  indicate  that  the  time- 
operator  splitting  approach  with  the  DVODE  solver  represents  a  viable  method  for  reducing  chemical 
stiffness,  allowing  utilization  of  low-cost  explicit  and  semi-implicit  diagonal  schemes  for  the  chemical 
transport  fractional  step. 


3  Accomplishments/New  Findings 

3.1  GASP  Governing  Equations 

GASP  uses  a  cell-centered  finite  volume  approach  to  solve  the  integral  form  of  the  Reynolds- Averaged 
Navier-Stokes  equations: 

h  ill Q  ^ + i(F  - FJ  -AiS= III  ™dV-  « 

The  standard  fluid-dynamic  system  has  been  augmented  to  include  the  effects  of  generalized  chem¬ 
istry,  i.e.,  the  global  continuity  equation  has  been  replaced  with  N  species  continuity  equations  to 
allow  finite-rate  chemical  processes.  The  conserved-variable  vector,  Q,  the  inviscid  and  viscous  flux 
vectors,  F  and  F„,  and  the  source-term  vector,  W,  are  given  as: 
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Here,  T  is  the  stress  tensor. 

For  a  structured  curvilinear  coordinate  system  (£,p,  (),  the  cell-centered,  finite  volume  formula¬ 
tion  for  Eqn.  (1)  can  be  written  as: 

tf+R(Q)  =  0  (2) 

where 

R(Q)  =  <%(F  -  F„)  +  <5„(G  -  G„)  +  <5C( H  -  H„)  -  i  W .  (3) 


3.2  Generalized  Chemistry 

GASP  allows  a  general  representation  of  chemical  systems  with  N  species  and  J  reactions  given  by 

v'l,jX\  +  V* IjX. 2  +  ...  +  V  N,j Xff  v" \  jX i  +  v" 2,jX-2,  +  ...  +  u" N,jX jv 

j  =  1, 2, ...,  J ,  (4) 

where  v's  ■  and  v”3  are  the  stoichiometric  coefficients  of  species  Xs  in  the  jth  reaction.  The  rate  of 
production  of  species  s,  denoted  us,  is  given  as 


(5) 


where  kfj  and  kf,j  are  the  forward  and  backward  reaction  rates.  In  GASP,  the  forward  rates  are 
determined  from  the  Arrhenius  equation 


kf  =  cT71  e~e/kT  , 


(6) 


and  the  backward  rates  are  determined  from  the  ratio  of  the  forward  to  equilibrium  rates 


kb  — 


kc 


(7) 


The  equilibrium  rates  are  determined  from  either  the  Arrhenius  equation,  equilibrium  curve  fits,  or 
from  the  McBride  Curve  fits  and  the  minimization  of  Gibbs  free  energy. 


3.3  Explicit  Stability  Bound 

Explicit  methods  for  solving  flows  with  finite  rate  chemical  kinetics  are  not  practical  because  the 
time-step  limit  imposed  by  the  explicit  stability  bound  is  too  small  to  make  the  computation  feasible. 
To  examine  this  fact,  we  consider  the  linear  scalar  advection  equation  with  a  “chemical”  source  term: 


dt  a  dx 


and  the  associated  explicit,  upwind  (a  >  0)  difference  equation 


A  tun  VxUi  u? 

— T—  -  = • 

A  t  Ax  tc 


(8) 

(9) 


Here,  rc  is  the  characteristic  time  for  chemistry,  and  is  representative  of  the  time  needed  for  the 
chemistry  to  reach  equilibrium.  Stability  analysis  produces  two  constraints:  one  based  on  the  wave 
propagation  (the  CFL  condition)  and  one  based  on  the  production  of  state.  Specifically,  we  have 


At  <  Ax/a, 
At  <  2  rc. 


The  second  constraint  requires  that  the  maximum  stable  time  step  decrease  with  decreasing  chemical 
time  scale  (rc) .  This  situation  is  labeled  “chemical  stiffness” ,  and  corresponds  to  flow  regions  where 
the  chemistry  is  near  equilibrium,  i.e.,  tc  —>  0.  The  chemical  time  scale,  for  a  given  reaction,  can  be 
related  to  the  reaction  rate  k  through  the  following  general  relation 


(10) 


where  u  is  a  function  of  species  concentrations.  Therefore,  rapid  chemical  kinetics  i.e.,  k  — >  oo, 
causes  numerical  stiffness  because  the  numerical  time  step  must  approach  zero  in  order  to  maintain 
stability. 


3.3.1  Implicit  Schemes 

In  general,  implicit  algorithms  bypass  chemical  stiffness  by  re-scaling  the  chemical  time  steps  to 
accommodate  the  fastest  reactions.  The  Euler  implicit  scheme  as  applied  to  Eqn.  (2)  is  written  as: 
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AQ  =  -R(Q), 


(11) 


Here,  A,  B,  C,  represent  the  inviscid  flux  vector  Jacobians,  Av,  B„,  Cv,  are  the  viscous  flux 
vector  Jacobians,  and  dW/dQ  is  the  source  term  Jacobian.  Eqn.  (11)  represents  a  banded  block 
(N  +  4)  x  (N  +  4)  matrix  equation  for  AQ,  where  N  is  the  number  of  chemical  species.  The  band 
width  is  dependent  on  the  choice  of  spatial  discretization  and  the  grid  size.  While  implicit  schemes 


provide  the  fastest  convergence  rate,  they  are  very  computationally  demanding,  both  in  terms  of 
memory  requirements  and  operation  count.  For  large  numbers  of  species  these  schemes  become 
computationally  prohibitive.  Fig.  1  shows  the  CPU  time  and  memory  required  for  the  lower-upper, 
symmetric,  Gauss-Siedel  (LU-SGS)  implicit  scheme  compared  to  the  Euler  explicit  scheme  on  a 
20  x  20  x  20  test  grid  using  GASP.  As  an  example,  the  CPU  time/iteration  for  60  species  is  888  s 
for  the  LU-SGS  algorithm  while  the  corresponding  CPU  time/iteration  for  the  Euler  explicit  scheme 
is  5.5  s,  a  factor  of  161.4  times  smaller.  The  required  memory  is  15.04  times  smaller  for  the  Euler 
explicit  algorithm.  Obviously,  there  is  great  incentive  to  utilize  explicit  or  semi-implicit  algorithms 
if  the  severe  time-step  limitation  can  be  removed. 

3.4  Loosely-Coupled  Chemistry 

The  new  time-integration  strategy  is  based  on  a  loosely-coupled  fluid-dynamic/chemistry  approach. 
This  method  allows  solving  the  5  fluid-dynamic  equations  ( i.e .,  the  mixture  continuity  equation, 
3  momentum  equations,  mixture  energy  equation,  and  two-equation  turbulence  if  included)  using 
implicit  methods,  with  the  possibility  of  solving  the  N  species  continuity  equations  using  low-cost 
explicit  or  semi-implicit  diagonal  schemes.  In  the  loosely-coupled  approach,  we  partition  the  gov¬ 
erning  equations  into  two  equation  sets;  a  chemical  set  and  a  fluid-dynamic  set.  The  partitioned 
conservative  variable  vector  is  given  as  follows: 


Here  Cj  is  the  specie’s  mass  fraction,  (i.e.,  c\  =  Pi/p).  In  the  first  step  of  the  loosely-coupled  procedure, 
we  update  the  fluid-dynamic  system  assuming  constant  mass  fractions,  Cj.  At  the  end  of  this  step  the 
species  densities,  p/,  are  adjusted  to  reflect  the  changing  mixture  density  as  pi  —  p'  ct.  In  the  second 
step,  we  solve  the  chemistry  system  (using  the  updated  flow  variables)  for  the  individual  specie’s 
densities.  During  this  step  the  mixture  density  is  held  constant,  and  the  species  mass  fractions  are 
determined  as  c*  —  pc%j  J2(p°i)-  This  summation  enforces  the  constraint  that  =  1-  'With  this 
methodology,  the  fluid-dynamic  system  handles  the  global  mass  conservation  while  the  chemistry 
system  specifies  the  chemical  composition.  Using  the  loosely-coupled  approach,  the  majority  of  the 
viscous  terms  are  handled  implicitly,  with  the  exclusion  of  the  species  diffusion  terms.  This  approach 
is  expected  to  adequately  reduce  the  viscous  stability  constraint  for  most  problems.  Our  approach 
also  allows  for  implicit  boundary  conditions  for  the  fluid  dynamic  system. 

In  our  recent  work  we  have  implemented  the  basic  framework  and  tested  the  loosely-coupled  scheme 
in  GASP.  Test  cases  show  that  the  method  convergences  well  and  produces  results  consistent  with 
the  fully-coupled  scheme. 
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(a)  CPU  Time  per  Iteration 


(b)  Required  Memory 


Figure  1:  CPU  time  and  memory  requirements  vs.  the  number  of  chemical  species  for  the  LU-SGS 
implicit  and  Euler  explicit  schemes  on  a  20  x  20  x  20  test  grid. 


3.5  Time-Operator  Splitting 

Solution  of  the  chemical  system  uses  a  first-order,  time-operator  splitting  approach  that  allows  iso¬ 
lation  of  the  stiff  source  term  from  the  transport  term.  The  chemistry  system  is  written  in  terms  of 
the  transport  and  reaction  components 

J^t1  +  Rct(Q)  +  Rcfl(Q)  =  0  ’  (12) 

where  Rcr  is  the  chemical  transport  residual  and  Rcr  represents  the  source  term  due  to  chemical 
reaction 

R-ct(Q)  =  -  FCJ  +  <5r;(Gc  -  GCJ  +  <5c(Hc  -  HCJ  ,  (13) 

Rc/t(  Q)  =  -jW.  (14) 

The  terms  are  then  split,  creating  two  systems  of  equations.  Thus,  a  single  step  of  the  first  order 
splitting  method  advancing  the  solution  from  time  tn  to  tn+1  =  tn  +  At  amounts  to  an  application 
of  time  discretizations  applied  to  the  system.  The  first  step  advances  the  transport  system 

j  +  rct(Q)  =  0,  (15) 

from  the  initial  condition  Qn  for  a  time  At,  and  the  solution  is  denoted  Qn+1.  The  second  step 
advances  the  reaction  system 

jir+Rc'i(Q)  =  0'  (16) 

from  the  initial  condition  Qn+1  for  a  time  At  to  yield  the  final  value  Qn+1. 

3.5.1  Transport  Fractional  Step  -  Explicit  Methods 

The  chemical  transport  fractional  step  (Eqn.  (15))  is  currently  solved  using  the  simple  low-cost  Euler 
explicit  scheme.  As  applied  to  Eqn.  (2)  the  Euler  explicit  scheme  is  written  as 

Q"+1  =  Qn  -  J  At  R(Q) .  (17) 

Future  work  utilizing  semi-implicit  diagonal  schemes  may  be  used  to  improve  stability  for  the  trans¬ 
port  fractional  step. 

3.5.2  Chemistry  Fractional  Step  -  DVODE 

The  chemical  reaction  fractional  step  (Eqn.  (16))  has  no  spatial  dependence  and  thus  is  essentially  a 
system  of  ordinary  differential  equations  at  each  node,  requiring  no  boundary  conditions.  The  chem¬ 
ical  production  fractional  step  is  solved  using  the  variable  coefficient  ordinary  differential  equation 
solver  DVODE  at  each  grid  node.  Input  to  the  DVODE  solver  requires  an  estimate  of  whether  or  not 
the  problem  is  stiff.  If  the  problem  is  stiff  DVODE  utilizes  Jacobians  (either  numerically  generated 
or  user-prescribed)  otherwise  if  the  problem  is  non-stiff  DVODE  utilizes  efficient  explicit  procedures. 
In  our  implementation,  we  use  an  estimate  of  the  local  chemical  Damkohler,  Da,  number  to  indicate 
chemical  stiffness.  At  locations  in  the  flow  field  where  the  Damkohler  number  is  larger  than  some 
threshold  value,  Damax,  we  specify  that  the  problem  is  stiff.  This  method  allows  DVODE  to  utilize 
Jacobians  only  where  necessary,  and  should  provide  maximum  efficiency  for  large  numbers  of  species. 
Development  of  methods  to  compute  the  Damkoher  number  are  discussed  in  the  following  sections. 
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Figure  2:  Fluid  particle  traveling  across  a  cell. 


3.5.3  Damkohler  Number 


The  Damkohler  number  for  chemistry  is  defined  as  the  ratio  of  the  fluid-dynamic  to  chemical  time 
scales 

Da= .  (18) 

Tc 

A  large  Damkohler  number  (Da  >>  1)  means  that  the  chemical  time  scale  is  much  smaller  than 
the  fluid-dynamic  scale,  i.e.,  the  chemistry  is  tending  toward  equilibrium  with  the  flow.  Likewise,  a 
small  Damkohler  number  means  the  flow  is  nearly  frozen,  and  chemical  changes  would  be  lagging  far 
behind  the  fluid-dynamic  changes.  The  local  fluid-dynamic  time  scale,  Tfd,  is  taken  to  be  the  time 
required  for  a  fluid  particle  to  advect  across  a  cell.  For  the  cell  shown  in  Fig.  2  the  fluid-dynamic 
time  scale  is  Tfd  =  Ax/V,  and  we  have 


Da  = 


Ax/V 

oj/k 


Large  Damkohler  numbers  indicate  fast  chemical  reactions  an  chemical  stiffness. 


(19) 


Chemical  Time  Scales 

Damkohler  limiting  is  based  on  knowledge  of  the  magnitude  of  the  chemical  times  relative  to  the 
convective  time  of  the  flow  (or  physical  time  step  for  unsteady  calculations) .  In  this  section  we  develop 
expressions  that  permit  estimation  of  the  characteristic  chemical  time  scales  for  the  elementary 
reactions  that  compose  a  chemical  kinetics  model  (see  Turns  [3]).  Notice  that  all  chemical  time 
scales  can  be  written  in  the  general  form  given  by  Eqn.  (10). 

The  simplest  possible  reaction  mechanism  is  the  unimolecular  reaction  given  by 

A  products  .  (20) 

This  reaction  occurs  spontaneously  by  decomposition  of  a  single  molecule  that  has  previously  been 
activated  to  a  high  energy  level  by  collision  or  otherwise.  The  corresponding  rate  equation  for  this 
reaction  is 

=  -k[A]. 


dt 


(21) 


Eqn.  (21)  can  be  integrated  to  yield  the  following  result  for  the  temporal  decay  of  [A] 


M 

Wo 


=  exp (— fct) , 


(22) 


where  [A]o  is  the  initial  concentration  of  specie  A.  The  time  constant,  rc,  is  defined  as  the  time 
elapsed  for  [A]  to  fall  from  its  initial  value  to  1/e  times  its  initial  value.  Hence,  we  solve  Eqn.  (22) 
for  the  time  when  [vl]/[v4]o  =  1/e,  yielding  the  following  chemical  time  constant 


r,  =  -  • 


k 


(23) 


As  an  example  of  a  bimolecular  reaction  which  occurs  at  a  collision  between  two  molecules,  we 
consider  the  following  reaction 

A  +  B  products  (24) 

and  its  rate  equation 

d\A] 

-k[A][B\.  (25) 


dt 


We  substitute  the  defect  variable,  x  =  [A]  —  [A]o  —  [£]  —  [B] o,  into  the  rate  equation  and  integrate 
to  yield 

dX  =-kt  +  C.  (26) 


/ 


(x  +  [A]0)(a:  +  [H]0) 

The  left-hand  integral  can  be  determined  to  be 

dx 


/ 


(x  +  [A]o)(a;  +  [ B]q )  Wo  —  [£]( 


In 


[B] 

W 


(27) 


and  setting  this  result  equal  to  the  right-hand  side  of  Eqn.  (26),  we  can  solve  for  the  temporal  decay 
of  [A]/[B]  as  follows 

|4|  =  j^exp{(P]„-[B]„)*t>  ■  (28) 

This  expression  is  used  to  determine  the  time  constant  for  the  decay  of  the  specie  concentrations. 
Assigning  [A]o  to  be  the  smaller  of  the  two  concentrations,  we  solve  Eqn.  (28)  for  rc  by  solving  for 
the  time  when  [A]/[A]o  =  1/e  where  we  use  the  fact  that  [B]/[B}0  =  ([A]o/[jB]o)([A]/[A]o  -  1)  +  1. 
We  obtain  the  time  constant  for  the  bimolecular  reaction  as 


Tr.  = 


WWq/WM 

([A)0-[B}0)k 


where 

and 


C(x)  =  —  lnje  —  (e  —  l)x] 


£{[A]0/[S]o}  =  ln  [A]/[^10 


[B]/[B\ o 


In 


(e-1) 


Wo 

[B)o 


(29) 


(30) 


(31) 


\[A]/[A]0=l/e 

As  an  example  of  a  termolecular  reaction  which  occurs  at  a  collision  between  three  molecules,  we 
consider  the  following  reaction 

2  A  +  B  products  +  A 


(32) 


and  its  rate  equation 


-m2[B}. 


We  again  substitute  the  defect  variable,  x  =  [A]  —  [A]o  =  [-B]  -  [f?]0,  into  the  rate  equation  and 
integrate  to  yield 

f  dx  _  ,  . 


J  {x  +  [^4]o)2(rr  +  [B]0) 

In  this  case  the  temporal  decay  of  [A]/[.B]  becomes 

[A\  [A]0 _ rDl  v2  1-4  ,  / 


=  -kt  +  C. 


This  expression  is  used  to  determine  the  time  constant  for  the  decay  of  the  specie  concentrations. 
For  the  case  where  [A]o  is  the  smaller  of  the  two  concentrations,  we  solve  Eqn.  (35)  for  rc  by  solving 
for  the  time  when  [A]/[A]o  =  1/e  where  we  use  the  fact  that  [B\/[B] o  =  ([A]o/[J3]o)([.<4]/[A]o  —  1)  + 1. 
We  obtain  the  time  constant  for  the  termolecular  reaction  as 

_  C{[A]0/[B}0}  -  X{[A]0/[B}0}  .  /ocN 


({A\o  -  [B]0)2  k 


where 


H*)  =  (X  “  (e“  x) 


[A}/[A}o=\/e 


Table  1  lists  the  chemical  time  scales  for  several  common  elementary  reactions. 

A  General  Method  for  Computing  the  Chemical  Time  Constant 

As  demonstrated  above,  the  analytical  expression  for  the  chemical  time  constant  depends  on  the 
form  of  the  elementary  reaction.  Rather  than  devise  expressions  for  every  possible  reaction,  we  have 
developed  a  general  method  which  applies  to  every  case.  In  the  context  of  Eqn.  (4),  consider  the  jth 
chemical  reaction  in  a  general  chemical-kinetics  model 

u'i,jX l  +  ^  %jX2  +  •  •  •  +  v'njXn  #  v" \jX\  +  y" 2, jX2  +  •  ■  •  +  v" n,jXn-  ‘  (39) 


The  defect  law  for  each  specie  can  be  written  as  follows 

[Xi]  -  [Xrlo  [X2]  -  [X2]Q 


V2 ,j  U2,j 


[Xn]  -  [X/y]o 
VffJ-VNj 


We  can  compute  the  chemical  time  scale  numerically  by  solving  an  ordinary  differential  equation  for 
the  specie  with  the  smallest  negative  defect.  That  is,  we  assume  that  vij  =  v”  ■  —  v[  ■  <  0  and  that 
the  species  are  listed  in  order  of  increasing  defect 

[*l]o  ^  P^]o  [Xn]o  /41n 

M  Kjl  Wnj  r  1  j 

Under  these  conditions,  we  can  solve  the  following  ODE  for  the  decay  of  [Xi] 

^r-  =  "  u  [x2]-«  ■  ■  ■  pc*]-".*.  (42) 


t 


Reaction 


Example 


rc 


A  — >  products 
A  +  M  — >  products  +  M 

2  ^4  — 4  products 

2  A  products  +  yl 
2A  +  M-*  products  +  M 

3  A  —>  products 

3  A  -4  products  +  A 
A  +  B  —>  products 

^4  +  B  +  M  —4  products  +  M 
2  y4  +  B  — 4  products  +  A 
2  A  +  7?  — 4  products 


7|-4  2I 

7720  +  M  — >  77  +  0/7  +  M 

2  0-4  02 

2  02  -4  2  O  +  02 

2T/  +  M— >T/2  +  M 

3  7  -4  72  +  7* 

3  77  -4  772  +  77 

772  +  02  — 4  2  770 

iV  +  O  +  M— >7V0  +  M 

2  77  +  077  -4  7720  +  77 

2  T72  +  02  — 4  2  7720 


1 

k 

1 

[M]7 

e  —  1 
2[A]0”fc 

e  —  1 

Mo  & 

e  —  1 

2Mo[M]fc 

e2  —  1 

6M§fc 

e2  —  1 

4MP 

£(s) 

([A)0-[B}o)k  ■ 
£0*0 

(Mo-[B]0)[M]fe 

£(x)  —  ^(x) 
([A]o-[7?]o)2fc 

C(x')  —  T(x') 
({A}o  -  [B}0)i  k 


Table  1:  Chemical  time  constants  where  x  =  [A]o/[R]o  and  A  is  assigned  to  the  specie  with  the  smaller 
initial  concentration  where  interchangeable.  The  second  to  last  reaction  assumes  that  [A]o  <  [J5]q, 
and  the  last  reaction  uses  x'  =  [A\q/2[B\q  and  assumes  [+l]o  <  2[7?]o. 


The  defect  law  (Eqn.  (40))  can  be  used  to  write  all  the  specie  concentrations  if  Eqn.  (42)  in  terms 
of  [Xl].  A  similar  ODE  can  be  written  for  the  smallest  decreasing  specie  in  the  backward  reaction 
with  rate  coefficient  kbj. 


To  solve  Eqn.  (42)  for  the  chemical  time  scale  we  use  a  second  order  mid-point  method  to  determine 
the  time  when  the  specie  decays  by  a  ratio  of  1/e.  We  first  make  a  first-order  approximation  using 
the  Euler-explicit  method  for  an  initial  guess 


[Xi]*  =  [Xx]o  + 


'd[Xi] 


dt 


T*ch- 


(43) 


Solving  for  T*h  and  noting  that  [Xi]*/[Xi]o  =  1/e,  we  have 


^ ~ch 


(1/e-l). 


The  second-order  approximation  can  be  written  by  re-evaluating  the  right  hand  side  of  Eqn.  (42) 
using  the  mid-point  concentration  (he.,  [Xi]o  +  [Xi]q  (1/2 T*h)).  Our  chemical  time  scale  for  any 
elementary  reaction  can  then  be  written  as 


rch  — 


d[X  i] 
dt 


(1/e-l). 


(45) 


Fluid-Dynamic  Time  Scales 

The  fluid-dynamic  time  scale  is  determined  locally  for  each  cell  in  the  flow  field.  For  a  structured 
curvilinear  coordinate  system  (£,77,  £),  we  compute  a  characteristic  fluid-dynamic  time  for  each  co¬ 
ordinate  direction 


r/d?  2 
Tfd  2 


+  ]^i+l/2l 
+  l^i+1/21 


^+1/2  \ 
Vi,j,k  J 

Aj+l/A 
Vij,k  ) 


(46) 

(47) 


Tfd  £  ~  2  (  l‘“*— X/2I 


be— 1/2 


2  V  Vititk 

Here  u  is  the  contravariant  velocity  on  each  face  given  as 


+  l^fc+1/2 1 


A-k+ 1/2  \ 

Vij,k  )  ' 


(48) 


u  =  u-Vk  (k  =  £,r},C), 


(49) 


where  Vk  is  the  face  unit  normal.  The  fluid-dynamic  time  scale  is  taken  as  the  minimum  of  each  of 
the  directional  times 


rfd  =  min {TfdvTfdjl,Tfd<;}  . 


(50) 


3.6  Results 

3.6.1  Ethylene  Flow  over  Blunt  Body 

As  a  test  case,  we  model  the  supersonic  flow  of  ethylene  (C2H4)  around  a  blunt  body -using  the 
7-specie,  chemical-kinetics  model  of  Baurle  et  al.  [1],  A  31  x  26  mesh  is  generated  around  an  axi- 
symmetric  body  traveling  at  =  7.46.  The  free-stream  pressure  is  px,  =  42,  657  N/m2  and  the 


densities  which  correspond  to  an  equivalence  ratio  of  <j>  =  1.667  are 


Pc2h4  =  0.0404635  kg/m3 
po2  —  0.0830436  kg/m3 
Pn2  =  0.2734928  kg/m3. 

The  Baurle  reaction  mechanism  is  given  as 

C2Ha  +  02  ^  2  CO  +  2  H2 
2C0  +  02  ^  2  C02 
2H2  +  02  #  2  H20. 

The  resulting  flow  solution  is  shown  in  terms  of  temperature  contours  (Fig.  3(a))  and  H20  mass 
fraction  contours  (Fig.  3(b)).  As  the  flow  passes  through  the  bow  shock,  the  increased  temperature 
ignites  the  ethylene  mixture  and  results  in  a  reaction  front  indicated  by  the  H20  mass-fraction 
contours.  Fig.  3(c)  shows  contours  of  the  Damkohler  number  for  this  problem.  The  Damkohler 
number  number  is  largest  near  the  stagnation  point,  reaching  a  maximum  value  of  Da  =  2223.84. 
A  significant  portion  of  the  flow  field  has  Damkohler  numbers  exceeding  Da  —  100.  We  ran  two 
calculations;  a  fully  coupled  case,  and  a  loosely  coupled  case  with  operator  splitting  and  the  DVODE 
solver  for  the  chemistry  system.  Jacobians  were  utilized  in  the  DVODE  solver  for  Damkohler  numbers 
greater  than  Damax  =  25.  Utilizing  the  loosely-coupled  scheme  with  operator  splitting,  we  have  been 
able  to  obtain  solutions  to  the  ethylene  blunt  nose  problem  using  the  low-cost  Euler  explicit  scheme 
for  chemistry  transport.  Fig.  4  shows  good  agreement  for  the  water  mass  fraction  along  the  stagnation 
line  for  both  schemes.  Without  the  operator  splitting  approach,  the  Euler  explicit  scheme  applied 
to  the  complete  chemistry  system  diverged  for  all  tested  time  steps. 

3.6.2  Hydrogen- Air  Wedge 

In  a  second  test  case  we  consider  the  two-dimensional,  inviscid  simulation  of  Mach  M  =  5.189, 
premixed  hydrogen  air  flow  over  a  20  deg  wedge.  The  grid  size  is  101  x  61  and  we  utilize  the  7- 
species  H2-Ah  model  of  Drummond  [2],  In  the  present  study,  the  inflow  mixture  is  fuel-lean,  with 
an  equivalence  ratio  of  —  0.75.  The  free-stream  conditions  are 

pN2  =  0.2374175  kg/m3 
Po2  —  0.0720856  kg/m3 
Pn2  =  0.0068121  kgjrr? 

T  =  800  K 
V  -  3300  m/s 

p  —  93.825  kN/m2 
M  =  5.189. 

Again,  we  ran  two  calculations:  a  fully  coupled  case,  and  a  loosely  coupled  case  with  operator  splitting 
and  the  ODE  solver  for  the  chemistry  system.  The  calculation  was  run  with  Van  Leer  flux  vector 
splitting  and  third-order  spatial  accuracy.  Fig.  5  shows  contours  of  temperature,  H20  fractions,  and 
Damkohler  number  for  this  problem.  As  the  flow  passes  through  the  shock,  the  increased  temperature 
ignites  the  H2-Aiv  mixture.  The  flame  front  occurs  downstream  of  the  oblique  shock,  and  is- indicated 
by  the  the  increase  in  water  mass  fraction.  Ignition  occurs  a  significant  distance  (i.e.,  the  induction 


(a)  Temperature. 


(b)  H2O  Mass  fraction. 


(c)  Damkohler  number. 


Figure  3:  Contours  of  temperature,  H2O  mass  fraction,  and  Da  number  for  the  ethylene  blunt-body 
flow. 
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Figure  4:  H2O  mass  fractions  along  the  stagnation  line  for  the  ethylene  blunt-body  flow. 

distance)  downstream  of  the  shock  due  to  the  finite  characteristic  time  for  chemical  reaction  and 
high  flow  velocity.  This  causes  a  separation  between  the  shock  and  flame  front  called  the  induction 
zone. 

Both  solutions  are  essentially  converged  after  2000  iterations.  Fig.  6  shows  the  temperature  and 
H2O  mass-fraction  profiles  along  the  j  =  41  grid  line  for  both  schemes.  The  results  show  very 
good  agreement  between  the  loosely-coupled  fully-coupled  solutions.  Convergence  hangs  up  slightly 
for  the  loosely-coupled  scheme,  but  both  solutions  converge  nearly  three  orders  which  is  typical  of 
higher-order  solutions  (see  Fig.  7).  For  the  7-specie  case  the  fully-coupled  approach  required  7,508  s 
while  the  loosely-coupled  approach  required  4, 940  s,  a  speed  up  factor  of  1.52.  To  examine  potential 
savings  for  larger  chemical  systems,  we  scaled  the  chemistry  model  to  include  up  to  20  species  by 
adding  inert  species.  For  the  20-specie  case  the  fully-coupled  approach  required  40, 485  s  while  the 
loosely-coupled  approach  required  11,706  s,  a  speed  up  factor  of  3.46. 

3.6.3  2-D  COIL  Simulation 

As  a  final  test  case  we  have  run  a  2-D  COIL  simulation  using  the  loosely-coupled  scheme.  This 
simulation  involves  l2-He  injection  into  a  primary  stream  of  excited  singlet  oxygen.  The  problem 
was  run  with  Roe  flux  differencing,  third-order  accuracy  and  viscous  effects  including  molecular 
diffusion.  Mach  number  and  small  signal  gain  contours  for  this  case  are  shown  in  Fig.  8. 

3.7  Conclusions 

We  have  examined  an  number  of  steady-state  flows  and  compared  results  with  loosely-coupled/operator 
splitting  methods  to  fully-coupled  results.  The  result  compare  well.  The  primary  potential  of  this 
approach  is  in  the  area  of  inexpensive  explicit  schemes.  With  the  source-term  stiffness  removed, 
the  memory  and  extra  computational  expense  of  fully-coupled  implicit  schemes  can  be  avoided.  We 
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(a)  Temperature. 
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Figure  6:  Line  plots  of  temperature  and  H2O  mass  fractions  for  hydrogen-air  wedge  flow. 
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Figure  7:  Convergence  history  for  fully  and  loosely-coupled  schemes  for  the  hydrogen-air  wedge  flow. 


(a)  Mach  number. 


(b)  Small  signal  gain. 


Figure  8:  Contour  plots  of  temperature  and  small  signal  gain  for  2-D  COIL  simulation. 


have  been  able  to  obtain  solutions  with  finite-rate  chemistry  using  explicit  methods  for  the  chemical 
system.  While  these  solutions  require  more  iterations  than  implicit  calculations,  they  become  much 
more  efficient  as  the  number  of  species  increases.  Although  this  work  is  focused  at  reducing  the  com¬ 
putational  resources  for  chemical  laser  calculations,  these  methods  are  also  applicable  to  hydrocarbon 
fuels  which  can  involve  50  or  more  chemical  species.  Detailed  chemical  mechanisms  for  hydrocar¬ 
bon  fuels  are  continually  progressing  in  their  accuracy,  however,  their  direct  use  in  calculations  is 
prohibitive,  and  reduced  models  must  be  used. 
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